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Absorption in dipole-lattice models of dielectrics 

R. J. Churchillf] and T. G. Philbiifl 

Physics and Astronomy Department, University of Exeter, Stocker Road, Exeter EX4 4QL, United Kingdom 

We develop a classical microscopic model of a dielectric. The model features nonlinear interaction terms 
between polarizable dipoles and lattice vibrations. The lattice vibrations are found to act as a pseudo-reservoir, 
giving broadband absorption of electromagnetic radiation without the addition of damping terms in the dynam¬ 
ics. The effective permittivity is calculated using a perturbative iteration method and is found to have the form 
associated with real dielectrics. Spatial dispersion is naturally included in the model and we also calculate the 
wavevector dependence of the permittivity. 

PACS numbers: 42.25.Bs, 77.22.-d 


I. INTRODUCTION 

Macroscopic electromagnetism remains a central compo¬ 
nent in describing a large variety of interactions between light 
and matter. The electromagnetic response of macroscopic ma¬ 
terials is routinely captured by an electric permittivity, and the 
permittivity itself can often be fitted to a simple function of 
frequency that belies the complexity of the underlying micro¬ 
scopic physics. Typically the permittivity £{uj) of a dielectric 
is expected to take the form 1 — an/ 
corresponding to a series of broadened resonances with atten¬ 
dant absorption (H. A simple argument (H, based on a po¬ 
larizable particle with damping together with the assumption 
of a rarified material, leads to this formula for the permittiv¬ 
ity. It is however difficult to justify this result with a more 
realistic model El. In particular, when one requires the ab¬ 
sorption of light to emerge from the model without being put 
in by hand as a damping term, one faces some difficult calcu¬ 
lations O. A classic paper by Hopfield in 1958 O elucidated 
the connection between the interaction of light with a polariz¬ 
able material on the one hand, and an effective description of 
such interactions by a permittivity on the other hand. Hopfield 
proposed that a rather simple model of the light-matter inter¬ 
action would lead to the permittivity given above, without the 
need to impose damping by hand or assume a rarified mate¬ 
rial. Although Hopfield’s model is simple to state and express 
as a Lagrangian, the derivation of the resulting effective per¬ 
mittivity requires considerable effort and does not appear to 
have been carried out before. The aim of this paper is to ver¬ 
ify Hopfield’s conjecture and demonstrate that the formula for 
£{(jj) above can be derived from a simple, intuitive, classical 
model of a dielectric. 

The initial model considered by Hopfield in ID is that of 
a continuous polarization field linearly coupled to light. This 
does not give realistic results however because light is only 
absorbed at the resonant frequency of the polarization field, 
giving a permittivity that has a delta function in frequency as 
its imaginary part m The problem with this model is identi¬ 
fied by Hopfield as an insufficient density of final states for en- 
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ergy transitions from light into the matter degrees of freedom. 
Hopfield states that the real reason for absorption of light by 
materials is nonlinear interaction between the electric dipole 
moments induced by light and lattice vibrations, so his polar¬ 
ization field would have to be nonlinearly coupled to lattice 
vibrations in order to see realistic absorption behaviour. No 
calculations however are given in ||2l for such a model with 
nonlinear matter interactions. Subsequently Hopfield’s sug¬ 
gestion has been followed up in a quantum setting where exci- 
tons (the quantized polarization field) are nonlinearly coupled 
to lattice phonons, the latter providing an effective damping 
when light interacts with the former |[3-[7l. A full analysis in 
the quantum setting is very challenging whereas an appealing 
aspect of Hopfield’s proposal is that it can be carried out clas¬ 
sically, much as the textbook motivations 111 for the macro¬ 
scopic Maxwell equations are classical. Here we will clarify 
that the broad absorption of the form used in the standard per¬ 
mittivity expressions can be derived from elementary classical 
physics. 

In Hopfield’s scenario the lattice vibrations act as a reser¬ 
voir into which electromagnetic energy can dissipate. This 
subsequently led to another approach in which a phenomeno¬ 
logical reservoir consisting of a field of harmonic oscilla¬ 
tors at every frequency is linearly coupled to the polarization 
field il [3. This phenomenological reservoir (an uncountable 
continuum of harmonic oscillators) is meant to mimic the dis¬ 
sipative channel into lattice vibrations through nonlinear in¬ 
teractions proposed by Hopfield. The continuum reservoir, 
introduced by Huttner and Barnett fSl, allows dissipation and 
a realistic permittivity to be derived from a linear model but a 
detailed connection to microscopic physics is less clear. The 
continuum reservoir approach has proven to be a powerful tool 
in incorporating dissipation in a Lagrangian with linear cou¬ 
pling. It has been used to give a Lagrangian formulation of the 
macroscopic Maxwell equations for an arbitrary permittivity, 
without use of a polarization field CMa, and similarly to 
give a Lagrangian formulation of damped harmonic oscilla¬ 
tors generally El. This gives a classical and quantum de¬ 
scription of light in all absorbing and dispersive macroscopic 
media, where the required permittivity (and magnetic perme¬ 
ability) are incorporated through coupling functions in the La¬ 
grangian. Here we return to the considerations ^ that led to 
the reservoir approach, to further clarify the microscopic pro¬ 
cesses that allow the macroscopic behaviour to be so accu- 
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rately captured by a phenomenological reservoir. We seek to 
verify that linear coupling of light to electric dipoles which are 
in turn nonlinearly coupled to lattice vibrations leads to a per¬ 
mittivity well described by the simple textbook formula given 
above, a permittivity that in particular exhibits the broadband 
absorption characteristic of real materials. Our analysis will 
also capture the nonlocal response (spatial dispersion) of the 
material medium, in addition to its temporal dispersion, giv¬ 
ing insight into the wavevector dependence of the permittivity 
of dielectrics. 

Because Hopfield’s proposal views the material as a lattice 
of dipoles, the model considered here is applicable to non- 
metallic solid dielectrics with a regular lattice. The model 
is thus not appropriate for amorphous materials, liquids or 
gasses, though in practice the permittivity functions of all 
these states of matter often show similar features. 

Although our motivation is to elucidate the microscopic 
physics of natural materials, it is of course possible to con¬ 
struct experimentally a lattice of dipoles on macroscopic 
scales. Metamaterials based on arrays of dielectric nanoparti¬ 
cles provide an example ESlIISl. Such artificial materials are 
being actively investigated as low-loss alternatives to struc¬ 
tures with metallic components 12014^ . The internal struc¬ 
ture of the dielectric nanoparticles, however, is an important 
aspect of such materials and nothing corresponding to this in¬ 
ternal structure is included in our model. (An important aspect 
of such dielectric nanoparticles is that the size of the particles 
can be used to control a magnetic as well as an electric re¬ 
sponse 1201 .) The coupling between the nanoparticles in the 
array is important but usually this coupling is investigated for 
its effect in altering the resonance structure of the metamate¬ 
rial |[2T][22l. In contrast, the important effect of the dipole 
coupling in our model is to induce the lattice vibrations that 
are essential for obtaining broadband absorption in our calcu¬ 
lations. On the other hand, in the microwave regime a lattice 
of dipoles that can vibrate is feasible in principle and may be 
useful to explore the absorption and non-local response cap¬ 
tured by our model. 

The outline of the paper is as follows. In Sec.|n|we con¬ 
sider a simple ID microscopic linear model and derive an ex¬ 
pression for the permittivity. We discuss the failings of such 
a model and the need for a nonlinear interaction term as a 
pseudo-reservoir. This is added to the model in Sec. ffl and 
the equations of motion are solved in Sec. IV using a pertur¬ 
bative iteration procedure. In Sec.|V|we introduce a graphical 
representation of the perturbative solution, which is used in 
Sec. VI to find the effective permittivity of the medium. Nu¬ 
merical calculations for the frequency and wavevector depen¬ 
dence of the effective permittivity are given in Sec. |VII| and 
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II. LINEAR MODEL OE A DIELECTRIC 


We initially consider a simple ID linear model that will 
prove inadequate to capture the absorption of light by materi¬ 
als. This will demonstrate why nonlinear interactions need to 
be incorporated into the model. 


Consider an infinite chain of polarizable dipoles positioned 
at Xn = na, where a is the lattice spacing. This model can 
also be considered as the limit of the large but finite chain, 
described by the same calculations with some small approxi¬ 
mations. Each dipole moment is modelled as a harmonic os¬ 
cillator with a Lagrangian: 
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where ujq is the resonant frequency of the harmonic oscillators 
and Tm describes the dipole-dipole interaction for a symmetric 
and translationally invariant system. This model of a dielectric 
can be viewed as a discrete version of the Hopfield model of 
a continuous polarization field O. The second part of our 
model is a ID scalar field (/)(x, t), representing a projection of 
the vector potential A, with the Lagrangian: 


L(f) 




{x,t) - , 


( 2 ) 


where c is the speed of light. The dipole moment is coupled to 
the time derivative of the scalar field, representing the electric 
field E = — A. An additional feature is a spatially dependent 
coupling term between the dipole moment Pn and the field 
near the lattice site The function a{x — x^) is used to 
account for the finite size of the dipoles, which are the “atoms” 
in our model. The interaction Lagrangian is given by: 



dx a(x — na)(j){x)^ 


(3) 


where do is the coupling strength and a{x) is taken to be a nor¬ 
malized Gaussian function, with the convenient feature that in 
the limit cr ^ 0 it becomes a Dirac delta function: 

o^{x) = —=e 2 . 2 . (4) 

ayziT 

At this point we make a spatial Lourier transformation. Lor 
the field, this is given by: 

1 

4>{x) = ^ dk4>{k)e^’^\ (5) 

V 27r J—OO 

where /c is a continuum over all values. Lor the infinite 
medium, the following expressions are used: 



( 6 ) 

(7) 


where q is now a continuous variable over the first Brillouin 
zone — ^ < q < ^. The delta functions in Q contain terms 
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where q is displaced by an integer number of the reciprocal 
lattice vector 27r/a. As we are primarily interested in initial 
fields with A ^ a, we only consider the j = 0 term at this 
point. Due to the real nature of the initial variables, the trans¬ 
formed variables obey (j){—k) = (t)^{k). The Lagrangian for 
this system now takes the form: 


^ 0 = 2 


ir 

2 A 


dk 


- k^(p{k)(p{-k) 


( 8 ) 


Lp- ^ 


\ j [p{q)p{-q) - i^l{q)p{q)p{-q)] , ( 9 ) 

a 

1 /“OO n dL 

<Pp = -do\J—J dk J dqp{q)a{-k)^{k)S{k +q), 


( 10 ) 


where wo(</) is the dipole dispersion equation in the absence 
of I/ 0 P and calculated from 

^ OO 

^ 0 (</) =^0 - o E 

m>l 

OO 

=‘^0 - E TmCos{qam), ( 11 ) 

m>l 


where we will only consider nearest neighbour coupling, with 
^m >2 = 0. The equations of motion of the Lagrangian are: 


'(j){k^t) -f {ck)‘^(l){k^t) = doc^ 



a{k)p{k,t), 


( 12 ) 


p{k,t) + u!Q{k)p{k,t) = -do-^ ^a{-k)4>{k,t). (13) 

Using the Fourier transformation: 

4,{t) = -= / (ko(j>{u)e 
V 27r j —OO 


(14) 


with the property (t){—k^ —uj) = (l)*{k,uj), the equations of 
motion become: 




(j){k^uj) = iuj^^V^oi{k)p{k^uj)^ (15) 

V CL 


[ci;o(^) — p{k^ uj) = —iuj^^\/^a{—k)<j)(k^uj). (16) 

V CL 

Solving gives: 

p{k,uj) =ph{k^uj) — iuj-^^V^CL{—k)Gp{k^uj)(l){k,uj), 

V Ci 

( 17 ) 

where ph is the homogenous solution of p satisfying the equa¬ 
tion of motion in the absence of coupling: 


[wo(fc) - Ph{k,uj) = 0, 


(18) 


and the retarded Green function Gp takes the form 
"a; 2 (fc)-(a; + i 0+)2 

=p _ - _I- 

^o(^) “ 

j'jr 

^-j^[5{u}-u}o{k))-5{u} + u}o{k))], (19) 


where the pole at oj^ = u}Q{k) has been moved into the lower- 
half complex plane by introducing the infinitesimal positive 
value 0+. This ensures that the solution is dependent on the 
field at previous times and satisfies the Kramers-Kronig rela¬ 
tions. Substituting ( p^ into ( p~5] ) gives 


-^e(fc, w) 


(p{k,uj) = iuj-^^'/^a{k)ph{k,uj), 
v CL 


where the relative permittivity 5 is given by: 


( 20 ) 
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s{k^uj) =1 H— (27r\a{k)\‘^) Gp{k^uj) 

^ (27r|Q((fe)p) 

uj^ik) — (cc + i 0 +)^ 


( 21 ) 


Equation ( p^ shows that for a given k, the imaginary part 
of the permittivity is given by a Dirac delta function at the cor¬ 
responding resonant frequency uJo{k), when the mode {k,uj) 
of the field cj) lies on the p dispersion relation. As uJo{k) is pe¬ 
riodic ink, a{k) ensures e ^ 1 as /c ^ oc for any non-Dirac 
delta function spatial coupling. As noted above, we only in¬ 
cluded the j = 0 term in 0 so the result is not the exact 
solution. Nevertheless the exact expression for the permittiv¬ 
ity also has an imaginary part that is a delta function. 

In reality, the complex permittivity near a resonant fre¬ 
quency has a finite imaginary component over a range of fre¬ 
quencies. This behaviour is usually modelled by treating each 
dipole as a damped harmonic oscillator (DHO), modifying Gp 
to the form: 


Gp{k, Lu) 


1 

ujQ{k) — ’ 


( 22 ) 


where 7 is a damping term. However, recovering the DHO 
equations of motion from a Lagrangian presents challenges. 
The oscillator can be coupled to either a discrete Il23ti26ll or a 
continuous llT7]| phenomenological reservoir. The oscillator- 
reservoir coupling must however take a very specific form in 
both cases if damping of the form seen in is to be recov¬ 
ered. In practice this damping term is usually added by hand 
to calculations. As this kind of damping leads to the standard 
permittivity expression stated at the beginning of this paper, 
we wish our model to produce it from physically motivated 
interactions, without extreme fine tuning of coupling terms in 
the Lagrangian. 

Hopfield faced the same problem of infinite absorption at 
a single resonant frequency in a similar dielectric model 0 , 
consisting of the electromagnetic field linearly coupled to a 
continuous harmonic-oscillator field (a polarization field, or 
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exciton field in quantum language). He identified the problem 
as due to the linear coupling, which allows only coupling be¬ 
tween single modes due to wavevector conservation. In terms 
of second-order perturbation theory, the photon exciton 
process does not have a density of final states, and so no 
real transitions occur. He then suggests that three-body (and 
higher order) interactions are responsible for absorption, fo¬ 
cusing on the exciton exciton' -\-phonon process. Unlike 
the linear case, there is an additional degree of freedom in in¬ 
teractions, with a continuum of coupled modes with the same 
total wavevector. This gives a continuum of final states with 
real transitions as energy absorbed from the electromagnetic 
field is stored in the indirectly coupled phonons. The coupling 
of the excitons to a pseudo-reservoir with the same wavevec¬ 
tor but a range of energies was part of the motivation behind 
the Huttner-Bamett model 0. Hopfield states O that such 
nonlinear interactions give rise to an effective damping term 
such as appears in (^), however no calculations are given. 
We now consider such interactions in a classical model. 


III. NONLINEAR MODEL OF A DIELECTRIC 


We modify the model of the previous section by expanding 
the behaviour of the medium into two degrees of freedom: the 
dipole moments of the “atoms” Pn and the physical displace¬ 
ment of the atoms from the equilibrium positions Xn given by 
Un (see Fig.[T]l. 

The atom displacement is not directly coupled to the field, 
instead being indirectly coupled through a nonlinear interac¬ 
tion with the dipole moments. This nonlinear coupling arises 
naturally from the dipole-dipole interaction: 


PnPn—m 


Lpp —Vo X] U , \ / , m3 

n=-cx3 m=-c»,7^0 I (^n-m +| 


c» oo ^ 

E E ^ 


PnPn—m 


n= — oo m= —00,7^0 


1 + 


^n — m 


3 • 


(23) 


Performing a small u expansion gives a series of Lagrangian 
terms with increasing powers of u. The first term is linear and 
is included in Lp with r^. The first nonlinear term, consider¬ 
ing only nearest neighbour |m| = 1 interactions, is 


L 


ppu 



a ) 


-PnPn +1 


(24) 


In addition, the atomic displacements are linearly coupled to 
each other by a Lagrangian term Lu : 


CX3 

-7;mo E 




oo 

m>l 


(25) 


where mo is the mass of the atom and Krn are the interaction 
terms. After a spatial Fourier transform: 


1 f ^ 

Lu = -jmo Idq [u{(l)H-(l) - , 

a 

(26) 

where the new resonant frequency is calculated in the same 
manner as uJo{q) in O: 


OO 

= kq — /^mcos(gam), (27) 

m>l 


and we again only consider nearest-neighbour interactions 
with hvrn >2 = 0. The nonlinear Lagrangian term takes 
the following form in Fourier space: 


T — ' 

J^ppu — r- - 

V27ra 


r a r a r a 

I I ^dqidq2dq3[fo{qi,q2,q3)piqi)p{q2)u{q3)5{qi+q2+q3)], 

Mqi,q2,q3) = i[sm{q2a) - sin((g2 + ga)®)], 


(28) 

(29) 


where 7 = 6 V 0 /a^ and we have omitted reciprocal lattice vector terms present due to 0 for notational brevity. After evaluating 
the Dirac delta functions, the equations of motion become: 


(j){k,uj) =(j)h{k,uj) (^V^a{k)^ G cf){k, uj)p{k, uj), 

p{k^uj) =ph{k^uj) — i^y^OL{—k^ Gp{k^uj)(j){k^uj) 


+ 


/ n \ r'^ 


fi{k,qi)u{qi,u!gjp{k - - WgJ 


i{q,u!) =Uh{q,u!)+j ( 2^^572 ) ^u{q,i^) J d{qia) J 


/ 2 {q, qi)p{qi, )p{q-qi,u!- uig^) 


(30) 

(31) 

(32) 
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^n-3 ^n-2 ^n-1 ^ ^n+1 ' ^n+2 ^n+3 


FIG. 1: Schematic of the nonlinear model of a dielectric. The dipole moments pn are linearly coupled to the scalar field cj) and nonlinearly 
coupled to the atomic displacement Un from the lattice site Xn — na. In addition, the atomic displacements are coupled to each other, with 
nearest-neighbour coupling displayed here. 


where the integration variables have been rescaled to dimen¬ 
sionless values and 

hik, qi) = -2i[sm{ka) - sin((fc - qi)a)], (33) 

f 2 iq, qi) = i[sin((g - qi)a) - sin(ga)]. (34) 

We have also introduced the retarded Green functions and 
Gu, given by 


G^{k,Lo) = 


(35) 

{ckY — {u! + i0+)^ ’ 

Gu{q,ojg) = 

1/mo 

(36) 

n^{q) - (w + i0+)2’ 


and the homogenous solutions of the fields (ph and satisfy 
the equations 

- [ckf] 4>h{k,uj) = 0, (37) 

[Wq - = 0. (38) 

Equations ( p^ and ([Tg for Gp and ph remain unchanged. 
Here we note that the nonlinear terms have introduced a 
pseudo-reservoir to the p equation of motion ( [ST] ), where a 
single initial mode is coupled to a continuum of modes with 
the same total k and uj. 


IV. PERTURBATIVE SOLUTION 

We now present a method of calculating an effective permit¬ 
tivity for this nonlinear system by deriving a wave equation 
for (j) similar to ( |2Q| ). This is done using an iteration proce¬ 
dure, treating terms with small nonlinear coupling coefficient 
7 as a perturbation of the linear model. We first consider ( [^ , 
as u is not directly coupled to p. In our model, the expres¬ 
sion for u{q,uj) can be immediately substituted into ( [ST] ) as 
it is expressed solely in terms of the homogenous solution Uh 
and p. For a system with different nonlinear coupling, the 
right-hand side (RHS) of ( [^ may contain additional terms 
of u{qPuj'). For example, further expansion of gives a 
ppuu term in the Lagrangian and a ppu term in In this 
case all RHS u terms in ( [3^ are repeatedly iterated. After 
n iterations, all RHS terms up to the power of 7 contain 
only Uh and p, while terms still involving u are of the order 
777+1 Qj. higher. Removing the remaining u terms leaves an 
expression for u{q^uj), accurate up order 7 ’^. 


The same situation is found upon substituting this expres¬ 
sion into ( [3T] ), with p{k, uj) in terms of homogenous solutions 
{uh^Ph)^ the field p and additional RHS terms of p{kPuj'). 
These terms are repeatedly iterated using the new equation of 
motion for p, to give an expression solely in terms of Uh, Ph 
and (p accurate up to order n in 7. 

A nonlinear wave equation for p is found upon substitution 
of p{k,uj) into ( [30| ). The RHS terms can be split into three 
groups: those containing only homogenous solutions (uh, Ph 
and ph), those linear in p and those nonlinear in p. The 
nonlinear p terms can be used to find an effective nonlinear 
permittivity; these terms can also be used to analyse the re¬ 
emission of frequency-converted p waves from an absorbed 
incident beam. These nonlinear terms start with higher-order 
powers of 7 compared to the linear terms and are dependent 
on higher powers of p. As a result, we can consider these 
terms to be negligible for “weak” fields. This leaves the terms 
linear in p. After substituting p{k,uj) into the linear p 
terms up to are: 


: 3 ' 


p{k^uj) (27r|(+(/c)|^) Gp{k^uj)p{k^uj) 

+ Z[(/)]+0[7"+i] (39) 


where Z[p] is a linear functional of p, containing an integra¬ 
tion J dk' f duo' over p{k',uj') terms. 

To find an effective permittivity of the medium from 1 
we perform a slightly different iteration procedure to the pre¬ 
vious two equations of motion. The linear p terms in Z[p] 
can be further split into two groups: those in the same mode 
pk^uo) as the other p terms in (39) and those in a different 
mode (k' ^ k^uo' ^ uo). The terms in the mode {k,uo) how¬ 
ever are a set of measure zero in an integration f dk' f duo' 
over all modes. This is a result of the continuous nature of 
q in •HJi for an infinite chain of atoms. The current model of 
an infinite chain must be treated as an approximation to the 
more realistic finite chain of N atoms. In the latter case the 
wavevector is a finely spaced set of N values. The integral 
over modes in the nonlinear process then becomes a discrete 
sum, where it is acceptable to separate a single term in the sum 
from the other terms. Thus, for the purposes of this iterative 
calculation it is necessary to treat the integration over modes 
as a sum over discrete values, whereas for numerical evalua¬ 
tion of final results the integral can be used without any appre¬ 
ciable error for a very long but finite chain. This is similar to 
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situations in quantum optics where a discretisation of modes 
renders some calculations easier, for example the treatment of 
thermal radiation fTfl . The RHS terms of 0 in ( [^ (including 
those from Z[(j)\) in the mode (/c, cc) are set aside, and terms in 
other modes are iterated using ( [^ . After n iterations we have 
a series of terms containing (j){k^uj) and powers of 7 up to 7 ^^, 
with leading order 7 ^ from the linear coupling; terms of or¬ 
der 7 ^^+^ and higher still contain (j) in different modes. The 
latter terms are dropped for an approximation to order 7 ^^. 
(Formally, if the iteration process is repeated indefinitely an 
expression with just (j) in the mode (/c, uj) will result.) When 
the iteration process is terminated and terms of order 7 ^+^ 
and higher are dropped, the resulting equation can be written 

[G(i){k^ (j){k^ uj) = iuj—^V^a{k)h{k^uj) 

V d 


+ {27r\a{k)\‘^) G'^{k,uj)(j){k,uj), 

G'Jk, Lo) = Gp{k, id) + 0(7), 


(40) 

(41) 


where h{k,uj) is now a collection of homogenous solution 
terms and Gp is now the leading order in a perturbation series 
giving the new function Gp. We identify G'p as an effective 
Green function describing the dressed dipoles p in the nonlin¬ 
ear medium. We can rewrite this equation in a form similar to 

[G'Jk^uj)] ^ (j){k^uj) = iuj^^V^a{k)h{k^uj)^ (42) 

V d 

where G'^ is the modified Green functions of 0, which can be 
used to find an effective linear permittivity 5eff: 

^ (Cfc)2-w2^eff(fc,w)’ 




s,«{k,co) = 1 + ^ (27r|a(fc)7 G'p{k,co). 


(44) 


V. DIAGRAMS 


by a pp term that can be iterated further. As a result all in¬ 
termediate lines, as shown in Fig. give a factor of the cor¬ 
responding Green function, while homogenous solutions that 
cannot be iterated further are represented as terminated lines. 


G^{k,uj) 


Gp (/c, cj) 


Gu{q,uj) 


-K -X 

(ph{k,uj) Ph{k,uj) Uh{q,uj) 


FIG. 2: Diagram representation of Green functions and homogenous 
solutions in the iteration process. 

The lines in a diagram may be connected with a limited 
number of allowed vertices, determined by the type of cou¬ 
pling in the Lagrangian. For example, the linear coupling 
in ( [ST] ) gives a two-line vertex, while the nonlinear coupling 
gives a three-line vertex. Each vertex in a diagram has an as¬ 
sociated prefactor from the equations of motion. The vertices 
and prefactors for the current model are given in Fig. At 
each vertex the total frequency of outgoing fields is equal to 
that of the ingoing field. Due to the periodicity of the mate¬ 
rial the wavevectors of p and u lie in the first Brillouin zone 
whereas the wavevector of (p has no such restriction. The to¬ 
tal wavevector of the outgoing fields at each vertex is equal 
to that of the ingoing field up to multiples of the reciprocal 
lattice vector. At a nonlinear vertex, all possible values of the 
undetermined qi and must be integrated over. The main 
difference to QFT is that we cannot use Wick contraction to 
close loops and remove additional pairs of u,por p terms. As 
a result, only tree diagrams are permitted. 

In summary: 

• Each intermediate line gives a factor of the correspond¬ 
ing Green function. 

• Each terminated line gives a factor of the corresponding 
homogenous solution. 

• Each vertex gives a factor of the corresponding cou¬ 
pling function from the equation of motion. 


The iteration procedure of the previous section rapidly be¬ 
comes notationally cumbersome. To simplify this process and 
the calculation of the effective permittivity, we express the it¬ 
erative procedure using diagrams. While Feynman diagrams 
were developed for quantum field theory (QFT) (see 1281 . 
for example), there is nothing inherently quantum about rep¬ 
resenting a perturbative solution to coupled field equations 
graphically. Feynman rules for diagrams can also be found 
when solving classical field equations perturbatively (291. 

The diagrams are to be read left to right. After each step 
in the iteration process, each field is represented as a line: 0, 
p and u are represented as wavy, straight and dashed lines re¬ 
spectively. Upon iteration using an equation of motion, a field 
is replaced by the homogenous solution plus the Green func¬ 
tion multiplied by a term involving another field. For example, 
in ( [ 3 ^ u{q,uj) is equal to Uh{q,uj) plus Guiq^dj) multiplied 


• Frequency uj is conserved at each vertex. The total 
wavevector is conserved at each vertex up to multiples 
of the reciprocal lattice vector, with the restriction that 
the wavevectors of p and u lie in the first Brillouin zone. 

• An integral is performed over each undetermined fre¬ 
quency and wavevector variable. 

• Only tree diagrams are permitted. 

This diagrammatic representation gives an intuitive way of 
finding the modified Green functions in a coupled system, by 
performing a summation over all diagrams that start and end 
with the same field. This is very similar to the calculation of 
the self-energy in QFT. This method also simplifies the iden¬ 
tification of terms that can either be grouped together or are 
part of an infinite series. 
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4>{k+^j,(^) p{k,uj) 

[V^ai-k-^j)] 

p{k,Uj) —Or\^ (j)(^k+ ^j,Ul) 
p{q-qi,u}- WqJ 

\ 

■“(?!> <^91 ) 

p{q- qi,uj - LJq^) 


u{q,ui) 



Using the rules described, we can evaluate this sum, with 
each line corresponding to a Green function and each vertex 
giving a factor of the coupling function: 


G^{k,Lu)= Gpik,uj) 


+ Gp{k^ cu) 


f; u;^^(2nla(k+^j)l^) 


X Gif,{k + oj) 


Gp(k,uj). 


(45) 


Dividing by GpGp gives: 

[G^piKoj)]-" =[Gp{k,aj)] 

oo 

- E 

j=-oo 


-1 


,2 ^0 


u:^^i2n\aik+fj)n 


xG^{k+fj,co) 


(46) 


FIG. 3: The allowed vertices and the corresponding coupling factors 
for the nonlinear model. 


As a simple example, we derive G^, the Green function 
of p dressed with 0, using only the linearly coupled model 
from Sec. [n| by performing a summation over all diagrams 
that start and end with p{k^uj). In the absence of coupling, 
we are left with only the first term of the bare Green function. 
Including the additional diagrams with intermediate (j) steps 
gives an infinite series that can be expressed via the Dyson 
equation |(28l: 


G(i){k + to) 

G^{k,uj) =Gp{k,uj) Gp{k,uj) Gp{k,uj) 

G(f>{k + ^j, uj) G(f,{k + u;) 

Or\^r\^r\jO 

^Gp{k,uj) Gp{k,uj) Gp{k,uj) 


G(f){k + Lo) 

= Gp{k,uj) +Gp{k,uj) Gd{k,uj) 


Reciprocal lattice vector scattering has been included explic¬ 
itly in this expression, as although the periodicity of the sys¬ 
tem restricts the initial wavevector of p to the first Brillouin 
zone, the intermediate (j) steps are not bound to this condi¬ 
tion. Exactly the same result can be found by substituting 
GD into Using the approximation a{x) = 6{x) and 

a{k) = 1 /the sum can be evaluated to give: 


G^p{k,Lu) = 




. f uja X 
sin( —) 


(47) 


9 .uja. X 

^ cos( — )—cos{ka) 


where the pole prescription {oo + i0+) has again been omitted 
for notational simplicity. 

The new dispersion relation for the dressed p is shown in 
Fig.[^ The reciprocal lattice vector scattering has the effect of 
folding the dispersion relation of (j) back into the first Brillouin 
zone, giving additional branches as j ^ oc. The (i-function 
approximation of a{x) is accurate for small uo but may not 
hold at very large frequencies {ooajc ^ tt) where the disper¬ 
sion relation is repeatedly folded back into the first Brillouin 
zone and j becomes large. 


VI. EFFECTIVE PERMITTIVITY 


We now calculate the effective Green function Gp{k,uo) in 
®, which is the Green function of p dressed by the nonlin¬ 
ear interaction; it reduces to the bare Green function Gp{k, uo) 
when 7 = 0 so linear pcj) vertices only occur between non¬ 
linear vertices. The iteration procedure also ensures that in 
Gp{k,oo) only intermediate p modes that differ from the ingo¬ 
ing p mode (k^oo) can couple to intermediate p lines in the 
ingoing mode (k^uo) do not connect to (j) lines. The diagrams 
for Gp{k,uo) are thus those that start and end with the bare 
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FIG. 4: The dispersion relation for the dressed dipole modes p lin¬ 
early coupled to a scalar field 0 for a finite do (red) compared to the 
dispersion relations of the uncoupled p and 0 folded back into the 
first Brillouin zone (dashed black). 


Green function Gp{k,uj), have increasing number of nonlin¬ 
ear vertices, and increasing numbers of intermediate p lines in 
the ingoing mode {k,uj). Similar to the previous section, this 
can be written as an infinite series: 


G'{k,uj) -Gp{k,uj) 



Gp (A;, cj) 


Gp(k,uj 



4 f(/c,u;)V 


Gp(k, co) 



{F(k,u;)} 


Gp(k, co) “ 


= Gp(k,uj) + Gp(k,uj) 



{F(k,u;)): 


Gp(k,uj) 


w) = w) 1 1 + F(fc, w)Gp(fc, w) 

+ [F{k,oj)Gp{k,oj)] + ... 

= Gp{k^Lu) [lF{k^Lu)Gp{k^Lu)'\ (48) 

Here, the modified Green function G'^ is represented by a 
straight double line and F{k,uj) represents a sum of all di¬ 
agrams that start and end with p in the mode {k,uj), where the 
two outer vertices are nonlinear vertices and the outer p lines 
are removed. Dividing by GpGp gives: 

[G'p{k,oj)]~" = [Gp{k,io)]-^ - F{k,oj). (49) 


The complex effective linear permittivity in ( |44| ) is now given 
by 


5eff(/c, LU) — 1 


(dgcVa) (27r|Q;(fc)|2) 


ujQ{kY — uj‘^ — MeF(/c, uj) — iImF{k, uj) 

(50) 

Instead of an imaginary Dirac delta function as in ( [2T] ), we 
have in (50) a resonant peak (27r\a{k)\‘^) /aImF{k,uj) 
when ujQ{k) — uj‘^ — ^eF{k^uj) = 0. The complex function 
F{k^uj) can be expanded in terms of the number of nonlinear 
vertices in each diagram: 


F{k,Lo) = 7^F2(fc,w) +'y^F3{k,u}) + 7^F4(fc,w) + ..., 

(51) 

where the term Fn contains diagrams where p returns to the 
initial mode after n scattering processes involving nonlinear 
vertices and n — 1 intermediate steps. Figure [^contains all di¬ 
agrams that start and end with Gp and contain two nonlinear 
vertices and corresponding powers of 7. The dressed Green 
function Gp from •EH' has been used to sum over all possible 
diagrams where the intermediate p step couples to (j) and back 
any number of times due to the linear coupling term. The dia¬ 
gram rules can be used to find the corresponding function for 
each diagram, which will contain an integral over the possi¬ 
ble final modes {k — qi — q 2 ^ 00 — ujq^ — The F 2 term 
will be calculated from the diagrams in Fig.j^by isolating the 
diagrams that return to the initial mode (/c, cc). 


Uh{qi,^qi) 

¥ 



I 

- Gu(k - qi,UJ - UJq^) 

Gp(k,uj) 

F2u " - 

Gp{k -\qi - q2,LU - Uq^ - Uq^) 

Ph(q2,u;q2) 


FIG. 5: All diagrams starting and ending with Gp and containing 
two nonlinear vertices. Note that the dressed Green function Gp in 
( [47] ) is used due to the linear coupling term. Each diagram will give a 
contribution to F 2 (k^ uj)\ we will denote these contributions by F 2 p 
(top) and F 2 u (bottom), based on the intermediate step. 


We now consider the first diagram in Fig. Using the 
Feynman rules, the corresponding function is 
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Gp(A:,a;)|y J d{qia) 7 ( 2 ^ 5 / 2 ) ( 9 i, WgJ G^(A: - gi, w - WgJ 

X J d{'^) J d{q2a) 7( ^^^5/2 )/i(^~gi’92) (92, WgJ Gp(A: - gi - ^ 2 , w - 0 ),^ - (52) 

where we have integrated over all possible final modes. This expression contains the homogenous solution Uh, which must 
satisfy ([38|). This can be expressed as a delta function in frequency: 


Uh {q,i^q)c 
2Tra^^^ 


= Ai „)5 + A -(-„)6 


+ 


ft{q)c 


(53) 


where we have included the preceding factor from the nonlinear vertex to simplify calculations and we have expressed the delta 
functions in dimensionless variables. The product Uh {qi^ojq^)uh (^ 2 , ^< 72 ) gives four terms, however we only consider 

those where ujq^ = —ujq-y , as only these terms allow a return to the initial mode 00 ) and thus contribute to F 2 : 


'Uh(qi,OJqJc' 

'Uh(q2,0JqJc' 

27ra^G 

27ra^/^ 


= A{qia)A^{-q2a)S 
FA'^{-qia)A{q2a)S 


Substituting the first two terms of •HI into ( [ 5 ^ and evaluating the delta functions gives 


LUq^a Q{qi)a\ (cUq^a ^{q 2 )ci 


^ ^{qi) a\ ^ ( idq^a n{q2)a 


(54) 






d {qia) fi{k, qi)A {qia) Gp{k - qi,uj - fi(gi)) 


/ d {q2a) fi{k - qi,q2)A* {-q2a) Gp(k - qi - q2,L0 - ^l{qi) + ^{q2)) 

J —TV 


(55) 


plus another term with A{qa) A*{—qa) and Q{q) —Q{q). At this point we separate ( [55] ) into terms that contribute to 
F 2 by returning to the initial {k,uj) mode and those that do not. By considering the integrals as a sum corresponding to a long 
but finite chain of atoms, we pick out the q 2 = —qi contribution in the second integral. The expression ( [55] ) then reduces to 
Gp{k,uj) ["f‘^F 2 p {k^uj)] Gp{k^uj)^ where F 2 p{k, uj) is the desired contribution io F 2 {k^uj) and is given by 

F 2 p{k,u)= j d{qia)\fi{k,qi)f\A{qia)fGp{k-qi,u:-Q.{qi)) 

J —TV 

d {qia) \fi{k, qi)f \A {-qia)f G‘i{k - qi,uj + fi(gi)). 
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(56) 


The integrals in (56) may contain poles of the dressed Green function G^, which coincide with the dispersion relation in Fig.|^ 
To evaluate ( [5^ we must insert the (uj + i0+) pole prescription in the G^ expression •ED’ which has the effect of shifting the 
pole into the lower-half complex u;-plane. We perform a change in the integration variables from wavevector q to frequency 
r^(g'), where it is easier to evaluate the poles: 


f d{qa)= [ Pi^)d{^), 

J —TV J O 




d(^] 


d(qa) 


(57) 


where denotes an integral over the finite ranges — to — and to f^niax<^/c using the expression from 

The expression ( [56] ) for F 2 p now takes the form 

F 2 p{k,Lo)= [ d{^)\Mk,Q{n,))f\A{^)\%{^)G^p{k-Q{n,),co-n,) 

Jn 


+ [ d{^)\Mk,-Q{Q^))f\A{^)fp{^)G^p{k + Q{Q,),oj-n,), 

Jn 


(58) 


r 


where Q{fl) is the inverse function of f^(g'), with the proper- ties Q{fl) = Q(—Q) and Q(Q) > 0. 
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The functional form of the lattice amplitude \A(ytia/c) \ in 
( [5^ must be specified. It is natural to take the homogeneous 
solution for the lattice as a thermal state. The average ampli¬ 
tude of a classical harmonic oscillator in thermal equilibrium 
is inversely proportional to its frequency so we take 

(59) 


As in the linear case ( p^ , the pole prescription of the re¬ 
tarded Green function G^n (58) can be used to split the in¬ 
tegrals into real principal value integrals plus imaginary terms 
associated with the poles. The latter terms correspond to the 
dressed dipole mode lying on the dispersion relation of Fig.|^ 
and can be found analytically in terms of residues; they give 
the imaginary part of ([58]) as 


where Aq is dimensionless. 


ImF2p{k^uj) =iiT ^Res 

n 

+i7r ^Res 

m 


|/i(fc, Q(Oi))| 2 |A f p Gl{k - Q(ni),a; - l^i), 

|/l(fc, IA (-^) 1% (^) G^ik + fii), , 


(60) 


where f^poie,n is the pole in the fti integration (the range 
of which is described after (57)) and Res[f{z),zo] denotes 
the residue of f{z) at zq. The principal value integrals that 
give the real part of ( |58| ) must be calculated numerically. An 
additional check on this numerical calculation can be made 
by comparing with the result obtained by using the Kramers- 
Kronig relations on the imaginary part (60) of F 2 p. 


In addition to F 2 p, there is another contribution to F 2 , asso¬ 
ciated with the second diagram in Fig.[^ and which we label 
F 2 u- This contributes very differently, however, for the fol¬ 
lowing reason. The linear coupling between p and (p means 
the poles of the intermediate Green function in F 2 p, which 
coincide with the dispersion relation in Fig. g occur for al¬ 
most all frequency arguments in Gp. This gives a nonzero 
imaginary part ( [60| ) of F 2 p{k^uj) for nearly every uj (recall 
that the residues in this equation are for the finite frequency 
range described after ([57])). In contrast, the poles in the in¬ 
termediate Green function Gu in F 2 U (see second diagram in 
Fig-IUf occur when The corresponding disper¬ 

sion relation runs over a smaller, finite range of frequencies 
determined by ( [27] ). This smaller range means that the first 
few leading order F terms with intermediate u steps will not 
have a pole in the integral over intermediate modes for a large 
range of initial frequencies uo. For example, in the region of 
interest near the dipole resonant frequency cc = cjq, F 2 u has 
no imaginary component for the choice of model parameters 
made in the next section. For this reason we only consider 
diagrams with intermediate dressed p steps in the numerical 
calculations that follow. 


The F 2 p expression (58) contains complex conjugate pairs 
of the homogenous solutions and the vertex coupling func¬ 
tion. The integration over intermediate states gives a construc¬ 
tively adding quantity as a result of this. Higher order terms in 
F do not necessarily have such complex conjugate pairs, in¬ 
stead containing a mixture of homogenous solutions and ver¬ 
tex functions at different frequencies and wavevectors. Upon 
integration, these can interfere destructively. In calculating the 


higher order terms, we only retain diagrams that give complex 
conjugate pairs, which give constructive interference and the 
dominant contributions to F. The next terms that satisfy this 
condition belong to the F 4 group and are calculated using the 
same process as F 2 p. The diagram and corresponding expres¬ 
sion for the leading order term F/^p is given in the appendix. 
Both F 2 p and F/^p will be used in the following numerical cal¬ 
culations, as they are the leading contributions to F. 


VII. NUMERICAL CALCULATIONS 


For numerical calculations, we consider a lattice spacing 
a ^ SA and use the approximation a{k) = l/\/^ for small 
initial k values within the first Brillouin zone. The resonant 
frequency ccq is taken to be in the visible region at 477.44 THz 
(1.975 eV), which corresponds to the dimensionless quantity 
ujQajc = 0.003. The coupling term rxiajcf = 1 x 10“^ 
in (U) is chosen so that uJo{k) is approximately constant in k. 
The lattice dispersion relation Q{q) covers a typical frequency 
range for a solid of approximately 0.1 10 THz The 

dimensionless terms Aq and ^(ajcf are chosen so that the 
F sum is convergent and perturbation theory is valid. For the 
purpose of our calculations Aq and are 0.0003 and 

1 X 10“^ respectively. Our general formulas are not specific 
to these values and are valid provided that the convergent F 
sum condition is satisfied. 

The most important constant is the linear coupling term do, 
which determines the size of the gap near cc = cjq in the 
dressed p dispersion relation, shown in Fig. [^(larger do gives 
a larger gap). If do is too large, the finite frequency integral in 
F2p over intermediate modes for an initial frequency uj ^ ujo 
will not include any poles of the intermediate Green function 
Gp or the corresponding imaginary residue terms. If do is 
too small, the integral includes poles from both the upper and 
lower branches of the dispersion relation in Fig. [^ with the 
possibility of the imaginary residues cancelling each other. In 
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both of these cases, one of the higher order terms such as F/^p 
will dominate the imaginary part of the F perturbation series. 
For the purposes of our calculations, we consider the inter¬ 
mediate case with d^^/a = 0.003 , where the integral over 
intermediate modes for an initial frequency uj ujq only in¬ 
cludes the lower branch of the dispersion relation and lmF 2 p 
is nonzero. 

We first consider /c = 0, where F 2 p dominates F for the 
chosen values of Aq and 7; in this k = 0 case the F^p contri¬ 
bution is not significant. Figure shows the real and imagi¬ 
nary parts of F 2 p near the dipole resonant frequency uq. The 
peaks in the imaginary part occur when the pole in the inter¬ 
mediate Green function Gp lies on the “fiat” part of the lower 
branch of the dispersion relation in Fig. The shape of the 
peaks is determined by /i(/c,cc) and A{Qa/c) from ( [3^ and 
( [59| ). The value of ImF 2 p( 0 , uj) between the two peaks near 
ujQ is small, but nonzero. As expected, the real and imaginary 
parts obey the Kramers-Kronig relations. 



0.00290 0.00295 0.00300 0.00305 0.00310 

codijc 

FIG. 6: RqF 2 p (blue) and ImF 2 p (red) rescaled to dimensionless 
variables in the frequency region near ujo{k) for k = 0. The real and 
imaginary parts are related by the Kramers-Kronig relations. 



^^>^(Hz) 


FIG. 7: The imaginary permittivity ( [M) (blue solid) and the DHO 
model fit ( [6^ (red dashed) when k = 0. The resonant peak has 
moved slightly from the position of the Dirac delta function of the 
linear model by ~ 1 kHz. The calculated permittivity is extremely 
well fitted by the DHO model. 


From ( [50| ) the imaginary part of the permittivity is 


Im5eff(/c,cc) 


(dgC^/a) [ImF{k^uj)] 

[wo(fc)2 - _ MeF(fc, Uj)f + [ImF{k, oj)f ’ 

(61) 


which has a resonant peak when a;o(fc) — = ReF{k,u). 

The small value of 7 and the shape of ReF{k,uj) mean the 
resonant frequency will only be shifted slightly from uJo{k\ 
The size of the peak is determined by ImF{k^uj). Figure]^ 
shows the resonant peak in Im^eff for /c = 0. The small values 
of 7 and Aq give the peak a small linewidth of ^ 60 Hz and a 
large maximum value of2.4xl0^^. The peak is only slightly 
shifted from the resonant value by ^ 1 kHz, corresponding to 
a2xl0“^^ fractional shift. The peaks in ImF 2 p in Fig. also 
give features either side of the central resonant peak in Fig.|^ 
These are smaller than the central peak by several orders of 
magnitude as the resonant condition cCq (/c) — = MeF(/c, 00 ) 

is not satisfied. The extremely sharp peak in Fig. [7] is due 
to our use of perturbation theory with a very small value of 
the nonlinear coupling parameter 7. A larger 7 would require 
more terms in the perturbation series to be evaluated, which 
involve more complicated intermediate scattering processes. 
More realistic results for the permittivity would require con¬ 
sideration of a very large number of intermediate processes, 
and permittivity values comparable to those measured in real 
dielectrics may be beyond the scope of perturbation theory 
in our model. However, we have shown that Hopfield’s pro¬ 
posal is correct: nonlinear interactions between the dipoles 
and lattice vibrations act as a pseudo-reservoir giving an ef¬ 
fective permittivity with finite imaginary part. The result¬ 
ing functional form of the permittivity is also in line with 
Hopfield’s statement O that such nonlinear material interac¬ 
tions should produce a permittivity agreeing with the standard 
damped harmonic oscillator (DHO) model 


£d{^) = 1 + 


Ad 


(62) 


Figure shows a fit of the imaginary part of the permittivity 
to the imaginary part of ( [62| ), where Ad,ujd and take the 
values Ad = 1.43 x 10^s“^, ujd = — 6240 and 

7 ^ = 202 s“^. It is readily apparent that our model gives an 
extremely good fit to the DHO formula ( [^ commonly used 
to describe real dielectrics. This close agreement is due to the 
fact that the imaginary term in the denominator of ( [^ ( 7 d<^) 
and nonlinear model (50) (ImF{k,uj)) do not change signif¬ 
icantly over the width of the peak. This imaginary term can 
thus be approximated as a constant, reducing both expressions 
to the Lorentzian function. The real part of the permittivity 
agrees with the DHO model to a similar degree, obeying the 
Kramers-Kronig relations. This behaviour is not specific to 
the nonlinear coupling used here and including more complex 
nonlinear coupling terms in the Lagrangian also gives a reso¬ 
nant peak that is an excellent fit to a DHO model, albeit with 
different parameters. 
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cd 

00 

"o 



ka 



The full /^-dependence SiS ka ^ tt would require calculating 
many terms and is not pursued further here. 

As k increases, the peak in the imaginary part of the permit¬ 
tivity (50), shown in Fig.|^for k = 0, decreases and broadens. 
This effect can be modelled by the DHO formula ( [62| ) by re¬ 
placing the parameters with a power series in k, for example 


^D{k) = CdDO + 0)^2^^ + Ldjj/^k^ + 


(63) 


where only even powers are present due to the symmetry of 
the system. For the numerical values used in this section, a 
k‘^ expansion in (63) provides a good fit to the calculated per¬ 
mittivity for /c up to /ca = 0.1, with further terms in (63) re¬ 
quired for higher k. A k‘^ expansion in (^), giving a pAerm 
in the denominator of ( [62| ), was proposed by Hopfield and 
Thomas ED based on different considerations. Information 
on the spatial dispersion of materials is limited in comparison 
to temporal dispersion, and the former is usually treated as of 
minor importance (321. Nevertheless an understanding of spa¬ 
tial dispersion is essential for accurate predictions in the nano¬ 
optics of small particles (3314351 and also for the prediction of 
Casimir and thermal forces on an isolated object (36l . Re¬ 
sults from our model may help clarify how spatial dispersion 
in dielectrics operates over a significant range of wavevectors. 


IX. CONCLUSIONS 


FIG. 8: Top: 'y‘^ImF 2 p{k,ujo{k)) (blue) and'yHmF 4 p{k,ujo{k)) 
(orange) as a function of k. Bottom: ImF{k,ujo{k)) as a function 
of /c, comparing the first (F 2 p) term (blue) to the sum of the F 2 p and 
Fap terms (red). F 2 p is the dominant term for small k within the light 
cone. At larger wavevectors, the residues in F 2 p begin to cancel and 
higher order terms become dominant in the F perturbation series. 


VIII. SPATIAL DISPERSION 

We now consider a nonzero wavevector in F{k^uj) to in¬ 
vestigate spatial dispersion in the medium, i.e. the wavevec¬ 
tor dependence of the permittivity (^). The first plot in 
Fig. shows the behaviour of the imaginary pats of the 
two leading terms F2p and F^p of F at the dipole resonant 
frequency uJo{k), as a function of k. The second plot in 
Fig. shows ImF{k^uJo{k)) with just the F 2 p contribution 
included and then also with the F^p contribution added. For 
small k, F2p{k,uJo{k)) gives the dominant contribution to 
ImF{k,uJo{k)). In this case the only processes in the dom¬ 
inant F 2 p term that give a residue in ( [^ are those from an 
integral containing poles from the lower branch of the dressed 
p dispersion relation. There are two such poles, with oppo¬ 
site signs of wavevector. As k increases, the residue term 
from the pole with the opposite sign of wavevector to k in¬ 
creases, while the residue from the pole with the same sign 
of wavevector as k decreases and changes sign. The overall 
value of F 2 p decreases and can become negative depending 
on the model parameters and coupling, as is seen in Fig. 
The higher order F^p term of F offsets the negative F2p term 
in Fig. with F/^p then dominating the expression for ImF. 


We have developed a simple classical model of a dielec¬ 
tric that features nonlinear interactions between polarizable 
“atoms” and lattice vibrations. Our motivation was to verify 
the main claims of Hopfield (21 regarding this model. Results 
such as those presented here give a better quantitative under¬ 
standing of the mechanism of light absorption in dielectrics, 
and also provide information on spatial dispersion. The lat¬ 
tice vibrations act as a pseudo-reservoir into which electro¬ 
magnetic energy is dissipated and the resulting permittivity is 


closely approximated by the standard textbook formula ( [62| ). 
As well as dissipation of incident radiation into the medium, 
our model will also describe re-emission of radiation out of 
the medium once the lattice is excited. The latter process is 
not captured by the effective permittivity and is contained in 
nonlinear (j) terms in our perturbation procedure that were not 
analysed here. 

We note that in our classical calculations the lattice must 
already be excited (we chose a thermally excited state as the 
coupling-independent part of the lattice solution) in order to 
perform its reservoir role. A quantization of the model would 
presumably also give broadband absorption of light at zero 
temperature due to the zero-point energy of the lattice. 

The reservoir role of the nonlinearly coupled lattice is cap¬ 
tured at a phenomenological level by linear coupling to a con¬ 
tinuum reservoir of harmonic oscillators at all frequencies O . 
The continuum reservoir, linearly coupled to the electromag¬ 
netic field, is in turn sufficient to give a Lagrangian formula¬ 
tion of the macroscopic Maxwell equations for arbitrary ma¬ 
terials obeying Kramer-Kronig relations (TOlfT^ . 

Spatial dispersion emerges naturally from our model. The 
wavevector dependence of the effective permittivity shows 
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agreement with simpler considerations ED for small k, but 
our model allows the calculation of higher-order contributions 
that are necessary for a full characterisation of the nonlocal re¬ 
sponse. 

The model explored here may also find application in mi¬ 
crowave metamaterials. Dipoles with a sharp resonance and 
very low internal loss that are arranged in a lattice that can 
vibrate could serve as a macroscopic system that is well de¬ 
scribed by our model. 
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APPENDIX: HIGHER ORDER F TERMS 


Uh{qi,UJq^) Uh{Q2,(^q2) 

X X 



FIG. 9: Diagram starting and ending with Gp and containing four 
nonlinear vertices, which gives a contribution to T 4 , which we denote 
F^p based on the intermediate steps. The additional restrictions on 
higher order F terms to give complex conjugate pairs of homogenous 
solutions and vertex coupling functions reduce the middle portion to 
the F 2 p diagram. 


The next term in the F sum after F 2 that satisfies the ad¬ 
ditional conditions in Sec. [Vl| (requiring complex conjugates 
of vertex functions and homogenous solutions) belongs to the 
F 4 group. While this contains many diagrams, we consider 
the term F^p calculated using the diagram in Pig.|^and named 
after the intermediate p steps. 


As before, the expression for the diagram in Pig. [^is found 
using the Peynman rules and the expression for F^p is found 
by only considering terms that return to the initial mode 
Due to the extra conditions we have imposed on the 
higher order F terms, we choose u{qs,ujq^) = u{—q 2 ^ “^< 72 ) 
and u{q 4 ,ujq^) = u{—qi^ ~^qi) to give complex conjugate 
pairs of the homogenous solutions and vertex functions. In 
this case the middle step of the diagram is now the same as 
The final expression for F^p can be reduced to 


F4p{k,co) = {^) I |A(fc,gKJ)|' |A {^)\G (- 5 ^) CO,,)]" 

X [F2p{k - Q{u!gJ,u - uigj] I 

+ ld{^)[\fi{k,-Q{cOg,))\^\A{-^)\"p{^) [G^^{k + Q{aogJ,ao-aogJ]" 
^ [F'2p{k ^ UJ ^qi)] r- 


The singularities in ( [64| ) are dealt with in the same manner as integral and a residue term, 
those in F 2 p, by splitting the integral into a principal value 


(64) 
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